program mainNxN
	
	implicit NONE
	
	integer(4), parameter				:: LX = 3
	integer(4), parameter				:: LD = LX*LX
	integer(4), parameter				:: MAXDIM = 2**LD
	
	integer(1)							:: repr(LD), rper(LD)
	integer(4)							:: upbound(LX), dnbound(LX)
	integer(4)							:: powers(LD)
	integer(4)							:: idd, iqq
	integer(4)							:: nd
	integer(1)							:: Hzm(MAXDIM,MAXDIM)
	integer(4)							:: Hex(MAXDIM)
	real(8)								:: Hm(MAXDIM,MAXDIM)
	
	integer(4)						:: sts
	
	real(8)							:: Jex, hx
	real(8)							:: evals(MAXDIM), E(MAXDIM), Z(MAXDIM), E2(MAXDIM)
	real(8)							:: emin
	
	real(8)							:: Zs, Es, Esqs
	
	integer(4)						:: lwork = 3*MAXDIM
	real(8), allocatable			:: wrk(:)
	
	integer(4), parameter			:: TMAX = 300
	integer(4)						:: eunit, tt
	real(8)							:: t
	real(8)							:: Eavr, Esqavr, Cavr
	
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	
	forall( nd = 1:LX)
		upbound(nd) = nd+1
		dnbound(nd) = nd-1
	end forall
	upbound(LX) = 1
	dnbound(1) = LX
		
	powers(1) = 1
	do nd = 2, LD
		powers(nd) = 2 * powers(nd-1)
	end do
	
	Hex = 0
	do idd = 1, MAXDIM
		call represent(idd,repr)
		call exdd(repr,Hex(idd))
	end do
	
	Hzm = 0
	do idd = 1, MAXDIM
		call represent(idd,repr)
		do nd = 1, LD
			rper = repr
			rper(nd) = 1 - repr(nd)
			call backrep(rper,iqq)
			Hzm(idd,iqq) = 1
		end do
	end do
	
!!!
	Jex = 1.0d0
	hx = 4.0d0
	Zs = 0.0d0
	Es = 0.0d0
	Esqs = 0.0d0

	Hm = - hx * Hzm
	
	forall( idd = 1:MAXDIM) Hm(idd,idd) = Hex(idd)

	allocate( wrk(lwork), stat = sts)
	call DSYEV( 'N', 'U', MAXDIM, Hm, MAXDIM, evals, wrk, lwork, sts)
	deallocate( wrk)


	emin = minval(evals)
	evals = evals - emin
	
	open( newunit = eunit, file='evals lx3h4.txt', status = 'replace', action = 'write')

	do tt = 1, 10 * TMAX
		
		t = 1.0d-2 * tt
	
		Z(:) = exp(- evals(:) / t)
		!forall(ip = 1:ND) Z(ip) = expSerL(- evals(ip) / t, 5)
		E(:) = evals(:) * Z (:)
		E2(:) = evals(:)**2 * Z (:)
		
		Zs = sum( Z(:))
		Es = sum( E(:))
		Esqs = sum( E2(:))
	
		Eavr = Es / Zs
		Esqavr = Esqs / Zs
		Cavr = ( Esqavr - Eavr**2) / t**2
		
		Eavr = Eavr + emin
	
		write (eunit, *) t, Eavr, Cavr
	
	end do
	
	continue
	
	close (eunit)

	stop
	
contains
	
subroutine represent( n, r)
	! Variables
	implicit NONE
	integer(4), intent(IN)			:: n
	integer(1), intent(OUT)			:: r(LD)
	!-----------------------------------------------------------------
	integer(4)				:: m, ip

	! Body of represent
	
	m = n-1
	do ip = 1, LD
		r(ip) = mod(m,2)
		m = m / 2
	end do

end subroutine represent

subroutine backrep( r, n)
	! Variables
	implicit NONE
	integer(1), intent(IN)			:: r(LD)
	integer(4), intent(OUT)			:: n
	!-----------------------------------------------------------------
	integer(4)				:: m, ip

	! Body of backrep
	n = sum( r(:) * powers(:)) + 1

end subroutine backrep
	
subroutine exdd( wf, mel)
	! Variables
	implicit NONE
	integer(1), intent(IN)			:: wf(LD)
	integer(4), intent(OUT)			:: mel
	!-----------------------------------------------------------------
	integer(4)			:: ip, iq
	integer(4)			:: x, y
	integer(4)			:: xt, yt

	! Body of exdd
	
	mel = 0
	
	do x = 1, LX
		do y = 1, LX
			ip = (x-1)*LX + y
			
			xt = upbound(x)
			yt = y
			iq = (xt-1)*LX + yt
			mel = mel + xor( wf(ip), wf(iq))
			
			xt = dnbound(x)
			yt = y
			iq = (xt-1)*LX + yt
			mel = mel + xor( wf(ip), wf(iq))
			
			xt = x
			yt = upbound(y)
			iq = (xt-1)*LX + yt
			mel = mel + xor( wf(ip), wf(iq))
			
			xt = x
			yt = dnbound(y)
			iq = (xt-1)*LX + yt
			mel = mel + xor( wf(ip), wf(iq))
			
		end do
	end do
	
	mel = mel - 2 * LD 

end subroutine exdd
	
end program mainNxN